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1. Introduction 

1.1 Imaging Atmospheric Cherenkov Telescope Arrays 

Imaging atmospheric Cherenkov telescope (lACT) arrays are the most sensitive instruments 
for detection of astrophysical gamma-ray emission at ~l-TeV energies [1]. They rely on the de¬ 
tection of Cherenkov light produced by particles in extensive air showers that were initiated by 
energetic astrophysical particles entering the atmosphere. The showers are imaged with multiple 
telescopes and their incoming directions and energies are reconstructed. Most of the astrophysical 
particles are hadrons, and while many of them are rejected based on the imaged shower’s proper¬ 
ties, some appear identical to photon-induced showers and are considered to be irreducible [2]. The 
dominant component in the count rate for most sources is this irreducible hadronic background; this 
is true even after reconstruction of the paiticle’s incoming direction, to an accuracy of about 0.1 
degrees. Thus, statistical tools for analysis must be employed and should be designed for maximal 
sensitivity in a background-dominated regime. 

The standard statistical method used by lACTs was introduced by Li&Ma in 1983. In their 
paper the following experiment is assumed [3]: an on-source observation is made for time To„ and 
the instrument is later shifted to observe a nearby background region (off-source) for time Toff- 
A test statistic is derived using maximum likelihood theory, given the ratio of observing times 
a = and the number of counts during the observations. Non and Noff- 

As lACTs evolved the tools to reconstruct the direction of events, the Li&Ma test statistic was 
generalized to the "wobble" analysis [4], where source (ON) and background (OFF) observations 
are made at the same time, but reconstructed to different spatial regions. The ring background 
model (RBM) is one of the standard tools for such analysis and has been successful in part due to 
its symmetrical nature that tends to cancel out systematic uncertainties in the background rate [4]. 

There are three current-generation lACT arrays, H.E.S.S. [5], MAGIC [6], and VERITAS 
[7]. CTA, a next-generation array with a sensitivity improvement of approximately an order of 
magnitude is expected to be in operation within a few years [8]. 

1.2 Gamma Ray Bursts 

Gamma-ray bursts were discovered accidentally in the 1960s by the Vela satellites [9]. While 
much progress has since been made on understanding their origin, the internal mechanism responsi¬ 
ble for their prompt radiation is still largely unknown [10]. The prompt emission is generally found 
in the keV-MeV range, but GeV emission has been observed by Eermi-EAT for several bursts [11]. 
Studies of Eermi-EAT light curves show an approximate power law decay of GRB fluxes [11]. 
Various models have attempted to explain the emission, some through synchi'otron self-Compton 
scenarios [12] which can involve electrons and positrons heated by a forward shock wave. Some 
of the proposed mechanisms imply substantial emission in the ~l-TeV energy range (for example, 
[13]). Gamma-ray bursts have been a coveted target for lACT observations, and are considered a 
primary target. None has been detected by any lACT array as of this date. 

The Ei&Ma likelihood method has been the standard approach for analysis of GRB obser¬ 
vations [14]. It uses the fact that the final count tally is Poisson distributed for both the source 
and background observations regardless of their particular time behaviour during the observation. 
Time varying sources can thus be inferred using this method, but any prior information on the time 
behaviour of these sources, as can be provided by other experimental observations or theoretical 
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predictions, cannot be included. We believe the inclusion of a priori temporal information can be 
an important tool in improving the sensitivity for GRB detection. 

Generalizations of the Li&Ma method have been introduced in the past [15]. They generally 
address the need to detect extended sources, and may need to rely on instrument response functions 
(IRFs) such as the point spread function or the energy reconstruction responses. To the best of our 
knowledge, no attempt has been made to explore a priori knowledge of the source light curve. 
We believe this can be achieved without the difficulties associated with the above methods, espe¬ 
cially since the time-stamping of events is much more accurate than the variability of any plausible 
gamma-ray source. 

The purpose of this paper is to discuss such a method. It is introduced as a natural generaliza¬ 
tion of the ring background model, and as such it is also resistant to systematic uncertainties and 
does not require detailed modelling of instrument response functions. 

Section 2 begins with a derivation of the Li&Ma test statistic (2.1), and then derives a simpli¬ 
fied form of a fime-dependenf fesf sfafisfic (2.2). The simplificafion relies on fhe assumption fhaf 
fhe hadronic background rafe is fime-independenf. Subsecfion (2.3) shortly describes a modifica¬ 
tion of the ring background model designed for adjusting to varying background rates. We believe 
this is the most robust approach for GRB detection. 

Section 3 describes results of an analysis of 8 promising VERITAS observations, along with 
an explanation for how those were selected. The analysis was performed with the method described 
in Section 2. 

2. Mathematical Derivation 

2.1 The Li&Ma likelihood ratio 

We will briefly derive fhe Li&Ma fest sfafisfic using average background and signal rales as 
free parameters inslead of counls (as used in fhe original paper). A more defailed discussion can be 
found in fhe original paper [3]. Our choice will serve as a smoolher fransifion fo fhe time dependenf 
fesf sfafisfic. 

In fhe likelihood model, fhe OFF counts are only due to an unknown background rate whereas 
the ON counts are explained by an unknown signal rate in addition to the same background rate. 

Defining fhe lime-averaged background and signal rales in relation lo fhe expected number of 
counls: bTo/f =< Noff >, {s-\-b)Ton =< Non >, the likelihood is given by: 


^ = P{No„\< Non >)P{Noff\ < Noff >) 

_ + b)Ton)^" {pTofff"" (2-1) 

To derive the null hypothesis likelihood, we simply set the signal rate to 0. In this case the 
average background rate is bo which completely accounts for all counts observed: 
b{) Poff — A V off ^ > bo Ton — ^ Non ^ • 

The null hypothesis likelihood is given by: 

,, {VoTonf- {VpTofff^" 

° Non'. Noff\ 
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We find the maximum likelihood values for the rates by maximizing the likelihood of both the 

No n +Noff _ _ Npff _ _ _ _ Npff 

Ton Tpff 


null and signal models: bg = 


Tpn + Tpff 

The likelihood ratio then simplifies into: 


b = 

Toff ’ 


■^0 

if “ (B + Sfonb^”" 


(2.3) 


Wilks’ theorem [16] allows us to describe the behaviour of the null likelihood ratio in the 

me of high counting statistics. If the null hypothes 
Gaussian variable with a standard deviation of 1 [16]. 


if 

regime of high counting statistics. If the null hypothesis is correct, \j —2log-^ is distributed as a 


2.2 Time-dependent signal, time-independent background 

To include arrival-time information, we divide Ton into an arbitrarily large number of equal 
bins N of time At, such that NAt = Ton- We will require the likelihood model to assign a probability 
for the number of counts within each bin independently. This will cause the likelihood to approach 
0 as —)■ oo because it will factor in the chance that the arrival times fall within specific bins, the 
number of which approaches infinity. This behaviour will cancel out in the likelihood ratio test, 
and the resulting test-statistic will converge nicely. 

In the limit of large N, each time bin will include either a single event, or no events at all. Each 
of the time bins is independently Poisson distributed with an expectation value approaching 0 as N 
increases. 

Let b denote the time-independent background rate, and s{t) denote the time-dependent signal 
rate. The background rate b is treated as an unknown to be optimized by maximum likelihood, 
and there may be similar unknowns within s{t), such as the amplitude or, a "shape parameter", etc. 
For use of Wilks’s theorem we must require nested models, and thus s{t) must have at least one 
such unknown, most simply the amplitude. If use of Wilks’s theorem is not possible, computer 
modelling can replace it, and the condition above can be relaxed. 

We denote the arrival times of signal (ON) photons as {ton} = {hd 2 --Tnp„)- The likelihood is a 
product of the Poisson probabilities for the count tally over all N time bins: 


{bToff)^<’ff ^-bTpff ^2.4) 

2At,...NAt) {0,1}! J Noff\ 

where [0,1] are chosen depending on whether there is an event in the f, bin. 

hm .5^ = i T\ {b + .(6))1 dtsb) 

V^,■€{/„„} / 

For the null hypothesis, we set s{t) = 0, and denote the background rate as bo, which will obey 
essentially the same likelihood ratio as the Li&Ma null hypothesis, with only a change of constants. 



ifo = At^^'-b^' 


{boToff)^‘’N bo(Ton+Tpff) 


( 2 . 6 ) 
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Thus ^0 is also maximized by: Zjq = giving: 




{boToff)^‘>ff (N,„+N,ff) 
Noff\ 


(2.7) 


The last equality follows from plugging in the exact form of Zjq into the exponential. The 
likelihood ratio is given by: 


±0 ^ _ b{T,„+Tff)+S^‘>’'dts{ty(Non+N,ff) 


( 2 . 8 ) 


This ratio can be further simplified by exploring the connection between s{t) and b. To do 
so, we must find the maximum of or equivalently of log .if. For the purpose of GRB detection 
we will leave only one free parameter, the amplitude, in the signal time profile, s{t) = df{t). 
This choice reflects the certainty of the flux decaying rapidly (usually as a power law), and fhe 
uncertainty about the amplitude of VHE emission. 


^logJ^" 

db 




1 

b+em 


{Ton + Toff) — 0 


(2.9) 


5logJ5f 

de 


I 

ti€{ton} 


fib) 

b + df{ti) 



( 2 . 10 ) 


An important identity is derived by noting that at the maximum of the likelihood function we 
can assert: b^^ + 6^^ = 0. 


b{Ton + Toff) + dt s{t) = Non +Noff (2-11) 

Jo 

The above relation leads to the cancellation of the exponential in the test statistic. It also 
allows a substitution of b = equation (2.9), which results in a polynomial 

equation of order Non + 1 for 0. It can be solved using a computer grid search, or any other suitable 
optimization algorithm to find a maximum likelihood 0. Due to the condition above, any grid 
search will only need to search over values of the amplitude, 0. 

We can now write a simplified form for fhe test statistic: 


C/> lJ^on-ri>o}} 

—9. = _0_ C2 12) 

In the case of a time-independent signal rate, this ratio is equivalent to the Li&Ma likelihood 
ratio, when written in terms of average rates. For variable light curves, ON counts with arrival 
times that match the expected profile of f{t) receive an essentially higher "weight", according to 
how bright the source was expected to be at the time. 
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2.3 A priori time-dependent background 

The background rate can also be modelled by an a priori time behaviour, for example by using 
a known dependence on the observation zenith angle, or using an independent part of the field of 
view for this purpose. This generalization of the previous section, based on the Ring Background 
Model [4], has been prepared and will be fully described in a future paper [17]. It is similar in nature 
to the technique described in Section 2.2 and is capable of handling a time-varying background rate. 

3. Analysis of VERITAS observations 

In the VERITAS observing program, gamma-ray bursts are considered the highest-priority 
sources. GCN alerts with a finer than 10° localization uncertainty radius are treated as a priority 
and prompt the observers to slew immediately to any such burst that is above a 20° elevation. 

VERITAS sensitivity is expected to be sufficient to detect strong bursts under reasonable as¬ 
sumptions about their high energy spectrum. Eor example, an extrapolation of EAT observations 
for GRB 130427A predicts an initial VERITAS flux of about 100 photons per second during the 
prompt phase [14], while, for comparison, the background rate for a region the size of the point- 
spread-function (PSE) is on the order of 1 count per minute. At the night of the burst, VERITAS 
was not operating due to a full moon; however, according to the estimates, a detection of GRB 
130427A would have been achieved in less than a second. Observations started the following 
night, and results yielded a small positive significance of 1.3 for the first night, and 1.1 for the 
second night [14]. 

As of the time this document has been written, VERITAS has observed 132 gamma-ray burst 
locations. Of those observations, we attempted to analyze the most promising candidates for de¬ 
tection based on a reasonable a priori measure. The most important factors in our ability to detect 
a burst are the observing delay, redshift (z), and elevation angle at the time of observing. The 
elevation angle is particularly important for high z bursts, since the instrument’s energy threshold 
rapidly increases with decreasing elevation angle, whereas absorption of gamma-rays due to the 
extragalactic background light (EBE) favours a very soft spectrum. We also required the burst to 
have been localized as well as our angular resolution or better (~ 0.1°). Bursts with unknown 
redshift were left out of the analysis, considering the high z typically associated with them. 

We calculated the following weight for each burst that met the conditions above: W = /tots, 

where T is the optical depth due to EBE absorption at the threshold energy, which depends on the 
elevation angle of observation 6, and Gfe is the observing delay. One can recognize this expres¬ 
sion to be approximately proportional to the incoming VHE flux at the initial time of observing, 
assuming the flux decays as 1/t, as is more or less typical for Eermi-EAT observed bursts [11] . 
The most promising observation was found to be that of GRB 150323A. We analyzed other bursts 
with weights within a factor of about 100 of its weight, to account for possible intrinsic variance in 
the VHE emission of the bursts. We found 7 other observations that met our criteria (see Table 1). 

The bursts were analysed using a standard VERITAS analysis package. Standard cuts opti¬ 
mized for soft spectrum sources were used in the analysis. The VHE flux was assumed to decay as 
a powerlaw with an index of — 1, beginning at the prompt phase of emission. The results for each 
of the 8 bursts is described in the Table 1. 

A form of a stacking analysis of the bursts was performed with what we will call the binomial 
test. It is performed as follows: Each burst is analyzed and its significance is found. A thresh¬ 
old significance is decided prior to analysis, and each burst can pass a test by having a greater 
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significance than this decided threshold, or otherwise fail the test. A p-value is calculated given 
the number of bursts that passed the test compared to what would be expected if their significance 
values followed a unit Gaussian distribution, as one would expect from the null hypothesis. 

The p-value is calculated with the following formula 
and then converted to Gaussian significance: 


r = (3.1) 

where I is the number of bursts passing the test, N 
is the total number of bursts in the analysis, and <I>(5) is 
defined here as the one-sided probability of obtaining a 
result with Gaussian significance greater than s: <I>(5) = 

The binomial test p-value may also be multiplied by 
the number of trials. In this analysis, we accounted for 
10 trials while using only 3 of them (see Table 2). The 
reason for this is that we are planning more trials on a 
wider selection of bursts without a known redshift. 


Name of Burst 

Significance 

GRB111225A 

-0.693 

GRB120422A 

2.473 

GRB130215A 

-0.600 

GRB130427A 

1.690 

GRB130604A 

-0.681 

GRB140622A 

-1.327 

GRB150120A 

-0.490 

GRB150323A 

-1.541 


Table 1: Time-dependent likelihood 
significance of 8 selected VERITAS 
bursts 


The threshold significance values were carefully selected using a computer script, to make 
sure that a whole number of bursts will be required to give an exact 5-sigma result. This resulted in 
fractional values for the thresholds, to a precision of 3 decimal digits. Table 2 summarizes our test 
selection (thresholds are calculated using a trials factor of 10). Table 2 describes the final p-value 
of all 3 binomial trials. 


Threshold 

significance 

Minimal number 

of bursts needed to 

achieve detection 

Number of bursts 
passing test 

p-value 

1.852 

6 out of 8 

1 out of 8 

0.23 

2.611 

4 out of 8 

0 out of 8 

1 

3.998 

2 out of 8 

0 out of 8 

1 


Table 2: Description and results of all binomial tests 


Summary 

We discussed a new technique for detecting time-dependent sources with lACTs; one that 
we believe can be especially useful for detection of gamma-ray bursts. We used this technique 
to perform a stacking analysis of 8 promising observations made by VERITAS. No significant 
detection is found. 
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